################################################################################
#  Replication materials for Graham & Orr (2020), "What Would Delegates Do?"
#  Study 1
#  This code generates: Figure 1, Figure A.1, A.2, and A.3
################################################################################


#------------------
#  Load packages
#------------------

library("mirtCAT")
library("ggplot2")
library("rstan")


#------------------
#  Define functions
#------------------

# Use 1D IRT framework to generate binary preferences
generate_pref <- function(difficulty, discrim, ideo){
  temp <- difficulty + discrim*ideo

  probvote <- 1/(1 + exp(-temp))

  if(probvote < 0) probvote <- 0
  if(probvote > 1) probvote <- 1
  
  return(sample(c(0, 1), prob = c(1 - probvote, probvote), size = 1))
}

# Use preferences and issue characteristics to estimate ideal points
issue_del_IRTideal <- function(Pref, model){
  issue_support <- apply(Pref, 2, mean, na.rm = TRUE)
  issue_support[issue_support == 0.5] <- NA 
  del_votes <- round(issue_support)
  return(unname(fscores(model, response.pattern = del_votes)[,"F1"]))
}

median_IRTideal <- function(Pref, model){
  observed_ideo <- fscores(model, response.pattern = Pref)[,"F1"]
  return(unname(quantile(observed_ideo, 0.5)))
}

# Graphics theme
theme_ideal = function(){
  list(theme_bw(),
       theme(
         panel.grid = element_blank(),
         axis.title = element_text(size = 9),
         legend.position = "bottom",
         plot.title = element_text(hjust = 0.5),
         strip.background = element_rect(fill = "gray10"),
         strip.text = element_text(color = "white")
       )
  )
}


#------------------
#  Figure 1, main text
#------------------

# Simulate ideal points varying ideological lean
set.seed(21519)

lean <- seq(from = -5, to = 5, length.out = 500)

sim_IRTdel <- rep(NA, length(lean))
sim_IRTmed <- rep(NA, length(lean))

K <- 50
J <- 2000

obs_discrim <- rnorm(n = K, mean = 1, sd = 0.25) 
obs_difficulty <- rnorm(n = K, mean = 0, sd = 0.5)

P <- matrix(NA, nrow = J, ncol = K)

for(i in 1:length(lean)){
  I <- rnorm(n = J, mean = lean[i], sd = 1)
  temp_model <- generate.mirt_object(data.frame(a1 = obs_discrim, 
                                                d = obs_difficulty),
                                     latent_means = lean[i],
                                     itemtype = '2PL')
  for(j in 1:J){
    for(k in 1:K){
      P[j, k] <- generate_pref(difficulty = obs_difficulty[k], 
                               discrim = obs_discrim[k],
                               ideo = I[j])
    }
  }
  
  sim_IRTdel[i] <- issue_del_IRTideal(P, temp_model)
  sim_IRTmed[i] <- median_IRTideal(P, temp_model)
}

#pdf("Figure1.pdf", width = 3.5, height = 3.5)
ggplot(data.frame(sim_IRTdel, sim_IRTmed), aes(x = sim_IRTmed, y = sim_IRTdel)) +
  geom_point(size = 0.25, alpha = 0.25) + 
  geom_abline(slope = 1, intercept = 0, color = "grey", linetype = 2) +
  labs(x = "Median ideal point", y = "Issue delegate ideal point") + 
  xlim(-5.15, 5.15) +
  ylim(-5.15, 5.15) +
  theme_ideal()
#dev.off()


#------------------
#  Figure A1 - A3, appendix
#------------------


# Difference from main text: Additive Score ideology estimates

# Define additional functions
issue_del_ideal <- function(Pref){
  issue_support <- apply(Pref, 2, mean, na.rm = TRUE)
  issue_support[issue_support == 0.5] <- NA
  del_votes <- round(issue_support)
  return(mean(del_votes, na.rm = TRUE))
}

median_ideal <- function(Pref){
  observed_ideo <- apply(Pref, 1, mean, na.rm = TRUE)
  return(unname(quantile(observed_ideo, 0.5)))
}

# Simulate new ideal points
sim_del <- rep(NA, length(lean))
sim_med <- rep(NA, length(lean))

for(i in 1:length(lean)){
  I <- rnorm(n = J, mean = lean[i], sd = 1)

  for(j in 1:J){
    for(k in 1:K){
      P[j, k] <- generate_pref(difficulty = obs_difficulty[k], 
                               discrim = obs_discrim[k],
                               ideo = I[j])
    }
  }
  
  sim_del[i] <- issue_del_ideal(P)
  sim_med[i] <- median_ideal(P)
}

#pdf("FigureA1.pdf", width = 3.5, height = 3.5)
ggplot(data.frame(sim_del, sim_med), aes(x = sim_med, y = sim_del)) +
  geom_point(size = 0.25, alpha = 0.25) + 
  geom_abline(slope = 1, intercept = 0, color = "grey", linetype = 2) +
  labs(x = "Median ideal point", y = "Issue delegate ideal point") + 
  theme_ideal()
#dev.off()


# Difference from main text: More controversial or more consensus issues

obs_difficulty_consensus <- rnorm(n = K, mean = 0, sd = 0.25)
obs_difficulty_controversial <- rnorm(n = K, mean = 0, sd = 0.75)

sim_IRTdel_consensus <- rep(NA, length(lean))
sim_IRTmed_consensus <- rep(NA, length(lean))

sim_IRTdel_controversial <- rep(NA, length(lean))
sim_IRTmed_controversial <- rep(NA, length(lean))

P_consensus <- matrix(NA, nrow = J, ncol = K)
P_controversial <- matrix(NA, nrow = J, ncol = K)


for(i in 1:length(lean)){
  I <- rnorm(n = J, mean = lean[i], sd = 1)
  temp_model_consensus <- generate.mirt_object(data.frame(a1 = obs_discrim, 
                                                d = obs_difficulty_consensus),
                                     latent_means = lean[i],
                                     itemtype = '2PL')
  
  temp_model_controversial <- generate.mirt_object(data.frame(a1 = obs_discrim, 
                                                d = obs_difficulty_controversial),
                                     latent_means = lean[i],
                                     itemtype = '2PL')
  for(j in 1:J){
    for(k in 1:K){
      P_consensus[j, k] <- generate_pref(difficulty = obs_difficulty_consensus[k], 
                               discrim = obs_discrim[k],
                               ideo = I[j])
      P_controversial[j, k] <- generate_pref(difficulty = obs_difficulty_controversial[k], 
                               discrim = obs_discrim[k],
                               ideo = I[j])
    }
  }
  
  sim_IRTdel_consensus[i] <- issue_del_IRTideal(P_consensus, temp_model_consensus)
  sim_IRTmed_consensus[i] <- median_IRTideal(P_consensus, temp_model_consensus)
  
  sim_IRTdel_controversial[i] <- issue_del_IRTideal(P_controversial, temp_model_controversial)
  sim_IRTmed_controversial[i] <- median_IRTideal(P_controversial, temp_model_controversial)
}

diff_robust_dat <- data.frame(del = c(sim_IRTdel_consensus, sim_IRTdel_controversial),
                              med = c(sim_IRTmed_consensus, sim_IRTmed_consensus),
                              robustness = c(rep("More Consensus Issues", 500),
                                             rep("More Controversial Issues", 500)))

#pdf("FigureA2.pdf", width = 7, height = 3.5)
ggplot(diff_robust_dat,
       aes(x = med, y = del)) +
  geom_point(size = 0.25, alpha = 0.25) + 
  facet_grid(.~ robustness) +
  geom_abline(slope = 1, intercept = 0, color = "grey", linetype = 2) +
  labs(x = "Median ideal point", y = "Issue delegate ideal point") + #title = "More Consensus Issues"
  xlim(-5.5, 5.5) +
  ylim(-5.5, 5.5) +
  theme_ideal()
#dev.off()


### Difference from main text: More or less ideological issues

obs_discrim_moreideo <- rnorm(n = K, mean = 1.5, sd = 0.25) 
obs_discrim_lessideo <- rnorm(n = K, mean = 0.5, sd = 0.25) 

sim_IRTdel_moreideo <- rep(NA, length(lean))
sim_IRTmed_moreideo <- rep(NA, length(lean))

sim_IRTdel_lessideo <- rep(NA, length(lean))
sim_IRTmed_lessideo <- rep(NA, length(lean))

P_moreideo <- matrix(NA, nrow = J, ncol = K)
P_lessideo <- matrix(NA, nrow = J, ncol = K)


for(i in 1:length(lean)){
  I <- rnorm(n = J, mean = lean[i], sd = 1)
  temp_model_moreideo <- generate.mirt_object(data.frame(a1 = obs_discrim_moreideo, 
                                                d = obs_difficulty),
                                     latent_means = lean[i],
                                     itemtype = '2PL')
  
  temp_model_lessideo <- generate.mirt_object(data.frame(a1 = obs_discrim_lessideo, 
                                                d = obs_difficulty),
                                     latent_means = lean[i],
                                     itemtype = '2PL')
  for(j in 1:J){
    for(k in 1:K){
      P_moreideo[j, k] <- generate_pref(difficulty = obs_difficulty[k], 
                               discrim = obs_discrim_moreideo[k],
                               ideo = I[j])
      P_lessideo[j, k] <- generate_pref(difficulty = obs_difficulty[k], 
                               discrim = obs_discrim_lessideo[k],
                               ideo = I[j])
    }
  }
  
  sim_IRTdel_moreideo[i] <- issue_del_IRTideal(P_moreideo, temp_model_moreideo)
  sim_IRTmed_moreideo[i] <- median_IRTideal(P_moreideo, temp_model_moreideo)
  
  sim_IRTdel_lessideo[i] <- issue_del_IRTideal(P_lessideo, temp_model_lessideo)
  sim_IRTmed_lessideo[i] <- median_IRTideal(P_lessideo, temp_model_lessideo)
}

discrim_robust_dat <- data.frame(del = c(sim_IRTdel_moreideo, sim_IRTdel_lessideo),
                                 med = c(sim_IRTmed_moreideo, sim_IRTmed_lessideo),
                                 robustness = c(rep("More Ideological Issues", 500),
                                                rep("Less Ideological Issues", 500)))

#pdf("FigureA3.pdf", width = 7, height = 3.5)
ggplot(discrim_robust_dat,
       aes(x = med, y = del)) +
  geom_point(size = 0.25, alpha = 0.25) + 
  facet_grid(.~ robustness) +
  geom_abline(slope = 1, intercept = 0, color = "grey", linetype = 2) +
  labs(x = "Median ideal point", y = "Issue delegate ideal point") + 
  xlim(-5.55, 5.55) +
  ylim(-5.55, 5.55) +
  theme_ideal()
#dev.off()


# Difference from main text: Probabilities generated with error

# Define additional functions
generate_pref_witherror <- function(difficulty, discrim, ideo){
  temp <- difficulty + discrim*ideo + rlogis(n = 1, location = 0, scale = 2)

  probvote <- 1/(1 + exp(-temp))

  if(probvote < 0) probvote <- 0
  if(probvote > 1) probvote <- 1
  
  return(sample(c(0, 1), prob = c(1 - probvote, probvote), size = 1))
}

# Simulate new ideal points
sim_IRTdel_error <- rep(NA, length(lean))
sim_IRTmed_error <- rep(NA, length(lean))

for(i in 1:length(lean)){
  I <- rnorm(n = J, mean = lean[i], sd = 1)
  temp_model <- generate.mirt_object(data.frame(a1 = obs_discrim, 
                                                d = obs_difficulty),
                                     latent_means = lean[i],
                                     itemtype = '2PL')
  for(j in 1:J){
    for(k in 1:K){
      P[j, k] <- generate_pref_witherror(difficulty = obs_difficulty[k], 
                                         discrim = obs_discrim[k],
                                         ideo = I[j])
    }
  }
  
  sim_IRTdel_error[i] <- issue_del_IRTideal(P, temp_model)
  sim_IRTmed_error[i] <- median_IRTideal(P, temp_model)
}

#pdf("sim_S_error.pdf", width = 3.5, height = 3.5)
ggplot(data.frame(sim_IRTdel_error, sim_IRTmed_error), 
       aes(x = sim_IRTmed_error, y = sim_IRTdel_error)) +
  geom_point(size = 0.25, alpha = 0.25) + 
  xlim(-5.15, 5.15) +
  ylim(-5.15, 5.15) +
  geom_abline(slope = 1, intercept = 0, color = "grey", linetype = 2) +
  labs(x = "Median ideal point", y = "Issue delegate ideal point") + 
  theme_ideal()
#dev.off()
